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Abstract 

A  continuous  adjoint  approach  for  obtaining  sensitivity  derivatives  on  unstmctured 
grids  is  developed  and  analyzed.  The  derivation  of  the  costate  equations  is  presented,  and 
a  second-order  accurate  discretization  method  is  described.  The  relationship  between  the 
continuous  formulation  and  a  discrete  formulation  is  explored  for  inviscid,  as  well  as  for 
viscous  flow.  Several  limitations  in  a  strict  adherence  to  the  continuous  approach  are  un¬ 
covered,  and  an  approach  that  circumvents  these  difficulties  is  presented.  The  issue  of  grid 
sensitivities,  which  do  not  arise  naturally  in  the  continuous  formulation,  is  investigated 
and  is  observed  to  be  of  importance  when  dealing  with  geometric  singularities.  A  method 
is  described  for  modifying  inviscid  and  viscous  meshes  during  the  design  cycle  to  accom¬ 
modate  changes  in  the  surface  shape.  The  accuracy  of  the  sensitivity  derivatives  is  estab¬ 
lished  by  comparing  with  finite-difference  gradients  and  several  design  examples  are  pre¬ 
sented. 
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Introduction 

Aerodynamic  design  optimization  has  been  an  important  area  of  research  for  many 
years.  Although  some  of  the  early  work  in  this  area  has  been  limited  in  applicability  be¬ 
cause  of  a  lack  of  computational  tools,  advances  in  computational  algorithms  and  com¬ 
puter  hardware  have  recently  fostered  intense  efforts  aimed  at  aerodynamic  and  multidis¬ 
ciplinary  optimization.  Among  the  methods  currently  used  are  gradient-based  optimizers 
in  which  a  specified  objective  function  is  minimized.  The  gradients  of  the  objective  func¬ 
tion  with  respect  to  the  design  variables  are  used  to  update  the  design  variables  in  order  to 
systematically  reduce  the  cost  function  to  arrive  at  a  local  minimum.  An  important  step  in 
this  process  is  the  determination  of  these  gradients,  which  are  also  referred  to  as  sensitiv¬ 
ity  derivatives. 

Several  techniques  have  been  investigated  for  evaluating  the  sensitivities  for  aerody¬ 
namic  applications.  A  description  of  these  techniques  can  be  found  in  Refs.  5, 15,  16,  and 
in  the  references  contained  therein.  Of  particular  interest  in  the  present  context  are  adjoint 
methods.  In  these  methods,  the  objective  function  is  augmented  with  the  flow  equations 
enforced  as  constraints  through  the  use  of  Lagrange  multipliers.  These  methods  are  partic¬ 
ularly  suited  to  aerodynamic  design  optimization  for  which  the  number  of  design  vari¬ 
ables  is  large  in  relation  to  the  number  of  aerodynamic  constraints  or  to  the  number  of  ob¬ 
jective  functions  in  a  multipoint  design.  This  is  because  the  derivatives  with  respect  to  all 
design  variables  for  each  objective  function  or  aerodynamic  constraint,  can  be  obtained 
with  a  computational  effort  roughly  equivalent  to  that  for  a  single  solution  of  the  flow 
equations. 

Adjoint  methods  can  generally  be  divided  into  discrete  and  continuous  adjoint  methods. 
In  the  discrete  adjoint  approach,  the  augmented  cost  function  is  discretized  before  varia¬ 
tions  are  taken.  For  the  continuous  adjoint  formulation,  the  process  is  reversed:  variations 
are  performed  first,  followed  by  the  discretization.  Note  that  the  operations  of  differentia¬ 
tion  and  discretization  do  not  commute  in  general.  Hence,  derivatives  obtained  by  using 
the  two  approaches  may  not  be  identical  and  would  differ  according  to  the  level  of  trunca¬ 
tion  error.  A  comparison  of  these  two  approaches  for  a  quasi-one-dimensional  problem  is 
given  in  Ref.  37. 

Much  of  the  pioneering  theoretical  work  in  adjoint  methodology  has  been  presented  in 
Refs.  19,  25,  29,  30,  and  31.  Although  optimality  conditions  for  aerodynamic  applications 
have  been  derived  from  a  continuous  approach  in  Refs.  3  and  6,  the  computer  implemen¬ 
tations  have  generally  followed  the  discrete  approach.  One  of  the  advantages  of  the  dis¬ 
crete  adjoint  approach  is  that,  because  the  equations  are  discretely  adjoint  to  the  flow 
equations,  the  derivatives  obtained  are  consistent  with  finite-difference  gradients  indepen¬ 
dent  of  the  mesh  size.  A  disadvantage  of  this  approach  is  that  it  requires  the  transpose  of 
the  matrix  that  represents  the  linearization  of  the  discrete  residual  with  respect  to  the  flow 
variables.  For  higher  order  accurate  schemes,  where  the  residual  has  a  complex  depen¬ 
dence  on  grid  points,  an  exact  implementation  of  this  approach  may  be  difficult  to  realize. 
For  this  reason,  previous  implementations  of  the  discrete  adjoint  approach,  such  as  those 
in  Refs.  6,  7,  26,  and  27,  have  used  a  discretization  of  the  adjoint  equations  that  is  consis¬ 
tent  with  a  first-order  accurate  discretization  of  the  flow  equations.  Second-order  accurate 
implementations  of  the  discrete  adjoint  approach  have  been  carried  out  on  structured  grids 
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in  Refs.  10  and  23.  On  unstructured  grids,  a  discrete  adjoint  approach  for  the  Euler  equa¬ 
tions  that  is  consistent  with  a  second-order  discretization  of  the  flow  equations  has  re¬ 
cently  been  implemented.*^ 

In  Ref.  19,  Jameson  developed  a  control  theory  framework  for  optimization  using  both 
the  full  potential  and  Euler  equations  for  compressible  flows.  Computational  results  based 
on  this  approach  were  first  presented  in  Ref.  20.  This  approach  has  been  further  developed 
and  implemented  for  both  two-  and  three-dimensional  applications.^*-^^  In  these  references, 
the  continuous  adjoint  approach  is  pursued  in  both  the  derivation  and  the  implementation 
on  structured  grids.  In  Refs.  32  and  33,  the  technique  has  been  applied  on  complex  config¬ 
urations  with  a  multiblock  algorithm. 

The  continuous  adjoint  approach  has  also  been  considered  by  lollo  et  al.*’  and  lollo  and 
Salas'*  for  both  one-dimensional  flow  and  two-dimensional  flows  over  simple  geometries. 
Kumvila,  Ta’asan,  and  Salas^'*  and  Ta’asan  and  Kuruvila**  have  investigated  an  efficient 
“one-shot”  approach  in  which  the  design  variables  are  updated  in  a  hierarchical  manner. 
Cabuk  and  Modi"  and  Cabuk  et  al.'^  have  also  used  an  adjoint  formulation  to  design  an 
optimal  diffuser  shape  using  the  incompressible  Navier-Stokes  equations. 

In  this  paper,  the  problem  of  aerodynamic  optimization  on  unstructured  grids  via  a  con¬ 
tinuous  adjoint  approach  is  developed  and  analyzed  for  inviscid  and  viscous  flows.  A  de¬ 
tailed  discretization  of  the  adjoint  equations  is  presented,  and  the  relationship  with  the  dis¬ 
crete  adjoint  approach  is  investigated.  The  accuracy  of  the  resulting  derivatives  is  assessed 
by  comparison  with  finite-difference  gradients.  In  addition,  a  mesh  movement  scheme  is 
presented  for  restructuring  the  grid  in  response  to  changes  in  the  surface  geometry.  The  re¬ 
sulting  methodology  is  then  used  to  design  several  airfoils  for  inviscid  compressible  flow, 
as  well  as  for  incompressible  laminar  flow. 

Adjoint  Variable  Approach  for  Sensitivity  Derivatives 

Considering  first  steady  inviscid  compressible  flow,  the  governing  equations  are  given 
by: 


|^F(Q)+|;G(Q)  =  0 

where  Q  is  the  set  of  dependent  variables  for  the  Euler  equations  (p,  pu,  pv,  E),  F  and 
G  represent  the  flux  vectors  of  mass,  momentum,  and  energy,  and  x  and  y  are  Cartesian 
coordinates. 

In  the  adjoint  approach  for  design  optimization,  a  cost  function  is  defined  and  aug¬ 
mented  with  the  flow  equations  as  constraints: 

I(Q,  D,  T)  =  I,(Q,  D)  +  R)dQ  =  I,(Q,  D)  +  Ir(Q,  D,  'R)  (2) 

where  R  represents  the  steady-state  flow  equations,  D  is  the  vector  of  design  variables, 
and  'R  are  the  Lagrange  multipliers  (also  referred  to  as  the  costate  or  adjoint  variables).  In 
Eq.  (2),  I  c(Q,  D)  represents  the  cost  that  is  to  be  minimized,  and  /^('R,  R)dQ  is  the 
inner  product  of  the  costate  variables  with  the  residual.  Examples  of  suitable  cost  func- 
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tions  include  drag  minimization  and  matching  a  specified  pressure  distribution,  for  which 
Ic(Q,D)  can  be  written  as 

Ij,(Q,  D)  =  ^^(cpkxcosa-cpkysina)ds  Drag  minimization  (3a) 


lx  2 

Ic(Q’ Specified  pressure  distribution  (3b) 


where  Cp  is  the  pressure  coefficient,  k^^  and  ky  are  x  and  y  components  of  a  unit  normal 
to  the  surface,  and  a  is  the  angle  of  attack.  The  cost  function  can  also  involve  field  inte¬ 
grals,  such  as  viscous  dissipation,  although  these  are  not  considered  in  this  paper.  It  is  as¬ 
sumed  that  the  cost  functions  are  differentiable  although  this  assumption  may  not  be  valid 
for  flows  with  shock  waves  or  other  singularities.  A  smoothing  procedure  as  suggested  in 
Refs.  19  and  22  may  be  employed  to  place  the  derivation  on  firmer  theoretical  ground. 
However,  in  numerical  implementations,  dissipation  typically  smears  discontinuities  over 
a  finite  number  of  mesh  points,  thus  mitigating  the  effects  of  non-differentiability.  There¬ 
fore,  smoothing  of  the  cost  function  is  not  performed  in  this  paper  with  no  apparent  conse¬ 
quences.  This  step  is  consistent  with  discrete  approaches  where  the  lack  of  differentiabil¬ 
ity  is  also  not  explicitly  taken  into  account. 

The  derivation  of  the  adjoint  equations  closely  follows  classical  techniques  from  calcu¬ 
lus  of  variations,  as  outlined  in  Ref.  36.  In  shape  optimization,  calculation  of  the  first  vari¬ 
ation  of  functionals,  such  as  those  in  Eqs.  (3a)  and  (3b),  requires  that  the  integral  on  the 
modified  surface  be  expressed  in  terms  of  quantities  on  the  original  surface.  For  example, 
cost  functions  such  as  drag  minimization  are  composed  of  terms  that  involve  products  of 
both  geometric  and  nongeometric  quantities: 

I,(Q,D)  =  f^(g(Q(D))k(D))ds  (4) 

Here,  g  is  an  arbitrary  function  of  the  flow  variables,  and  k  represents  either  k^^  or  ky . 
For  cost  functions  such  as  Eq.  (3b),  k  assumes  a  value  of  unity.  A  general  form  for  the 
first  variation  can  be  written  as 

6I,(Q,D)  =  ^^,(g„g^k„,^)ds  -f^(goidkoid)ds  (5) 

/ 

where  F  and  F  represent  the  old  and  the  new  surface  of  the  geometry,  respectively,  and 
the  subscripts  old  and  new  denote  quantities  on  these  surfaces.  Evaluation  of  these  inte¬ 
grals  is  addressed  after  a  discussion  on  obtaining  variations  of  Ij^ . 

The  method  for  obtaining  the  variations  of  the  volume  integral  in  Eq.  (2)  involving  the 
residual  R  follows  closely  that  of  Pironneau.^*  Denoting  this  volume  integral  as 
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(6) 


Ir(Q,D,'F)  =  J('F,R)dQ 
£2 

the  variation  that  properly  accounts  for  volume  changes,  as  well  as  for  changes  in  the  flow 
field,  is  given  by 

5Ir  =  -  J  +  B^0jdQ  +  J  Q^(  a\^  +  B\y)'Pds  (7) 

T  T 

where  A  and  B  are  the  transposes  of  the  inviscid  flux  Jacobian  matrices  and  the  sur¬ 
face  integral  is  over  the  solid  walls  as  well  as  the  far-field.  In  deriving  Eq.  (7),  it  is  tacitly 
assumed  that  the  fluxes  and  the  costate  variables  are  differentiable;  similar  assumptions 
have  been  discussed  earlier  regarding  the  cost  function.  The  variation  of  the  augmented 
cost  function  in  Eq.  (2)  is  formed  by  combining  Eq.  (7)  with  the  variations  in  I^CQ,  D) . 
Because  Q  is  arbitrary,  the  volume  integrals  present  in  the  variation  of  the  augmented 
cost  function  can  be  eliminated  by  requiring  that  'F  satisfy  the  following  adjoint  (costate) 
equation: 


-B  ^  =  0 
dy 


(8) 


The  surface  integral  in  Eq.  (7)  is  used  together  with  the  variations  in  the  cost  function 
Ij,(Q,  D)  to  determine  both  the  boundary  conditions  and  the  sensitivity  derivatives.  The 
boundary  conditions  for  *F  are  chosen  to  eliminate  the  terms  that  multiply  Q  on  the 
boundaries.  The  surface  integral  can  be  rewritten  as 

Jq\a\  +  B\y)Tdr  =  {(aV  Q)dr  (9) 

r  r 

— T  T  T 

where  A  =  A  +  B  ky .  In  the  far  field,  this  term  can  be  rewritten  by  using  a  locally 
one-dimensional  characteristic  decomposition  at  the  boundary  to  yield 

|(AVQ)dr  =  J('F,TAW)dr  (lO) 

r  r 

where  W  =  T~*Q,  is  the  matrix  of  left  eigenvectors  of  A,  and  A  are  the  corre¬ 
sponding  eigenvalues.  Boundary  conditions  for  the  costate  variables  in  the  far  field  are  ob¬ 
tained  using  characteristic-type  boundary  conditions  on  the  field  variables,  where  the 
propagation  of  information  is  based  on  the  signs  of  the  eigenvalues.  For  shape  optimiza¬ 
tion,  variations  in  W  associated  with  free-stream  quantities  are  zero,  so  that  the  corre¬ 
sponding  costate  variables  on  the  boundary  can  be  extrapolated  from  the  interior  of  the 
domain.  The  other  costate  variables  on  the  boundary  are  obtained  by  requiring  the  remain¬ 
ing  terms  in  Eq.  (10)  to  vanish.  When  Mach  number  or  angle  of  attack  are  design  vari- 
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ables,  variations  in  W  reflect  the  appropriate  changes  in  free-stream  conditions  and  are 
used  in  obtaining  the  derivatives  with  respect  to  these  variables. 

On  solid  walls,  the  boundary  condition  that  there  is  no  flow  normal  to  the  new  surface  is 
written  as 


y  +  Q2j(kx  +  kx)  (11) 

■^{,*53+5^  x+3^y  +  Q3j(k,+k,)  =  o 
Using  Eq.  (1 1),  the  surface  integral  in  Eq.  (9)  can  be  expressed  as 


J  ( A^'P,  Q)ds  =  |(Q  1  [(kx\|/2  +  ky\|/3)<t)]  + 
r  r 

Q2[(kxV2  +  ky\|/3)(l  -  Y)u]  + 

Q3[(kx¥2  +  ky\ir3)(l-y)v]  + 

Q4[(kxV2  +  kyV3)(Y  -  1)]  +  R(¥i  +  ¥2^  +  ¥3V  +  ¥4H)  )ds 


where  H  is  the  total  enthalpy  and 


R  =  -(Q2kx  +  Q3ky) 


^Q3V 
ydx  r 


(12a) 


(12b) 


In  order  to  compute  the  variation  in  Eq.  (5),  the  integrand  for  the  first  integral  is  ex¬ 
panded  as  follows: 


Snewkflew  =  (g  +  gxX  +  gyy  +  g)(knew  +  knew)  (13) 

In  Eq.  (13),  the  derivatives  g^  and  gy  account  for  spatial  changes  and  g  reflects  the  vari¬ 
ation  due  to  the  fact  that  the  solution  of  the  governing  equations  has  changed  in  response 
to  the  changing  surface.  Note  that  for  structured  grids,  which  employ  a  mapping  to  a  fixed 
computational  domain,  these  spatial  derivatives  do  not  arise  because  the  variations  in  the 
generalized  coordinates  are  zero.  However,  variations  in  the  mapping  function  need  to  be 
considered  which  naturally  provides  a  mechanism  to  account  for  grid  sensitivities  in  a 
continuous  framework.^^ 
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The  boundary  conditions  for  the  costate  variables  are  derived  by  combining  the  bound¬ 
ary  terms  from  the  variation  in  the  cost  function  with  those  from  Eq.  (12a)  and  then  elimi¬ 
nating  terms  that  involve  variations  in  Q .  Because  g  is  a  function  of  Q ,  g  is  given  by 


g  =  Y  |i-Qi 


(14) 


The  boundary  terms  that  multiply  Q  are  eliminated  by  requiring  that 


■  ■ 

(l-y)u 

+  k 

9g  3g  3g  9g 

(l-y)v 

.(y-1). 

3Q,’3Q2’3Q3’3Q4 

=  0 


(15) 


Note  that  the  column  vector  that  multiplies  (kx\K2  +  corresponds  to  the  derivatives 

of  pressure  with  respect  to  the  dependent  variables.  In  order  to  obtain  a  unique  boundary 
condition  for  (k^\j/2  +  kytitj) ,  the  second  column  must  be  a  scalar  multiple  of  the  first. 
Therefore,  g  can  only  be  a  function  of  pressure,  h(p) ,  which  yields  the  following  bound¬ 
ary  condition: 


k^\|r2  +  ky\|t3  +  k^h(p)  =  0  (16) 

Thus,  cost  functions  such  as  specification  of  a  velocity  distribution  or  minimization  of  sur¬ 
face  entropy  are  inadmissible,  except  in  special  cases  where  they  can  be  expressed  solely 
in  terms  of  pressure. 

As  an  example  of  an  allowable  cost  function,  consider  the  drag  coefficient  given  by 

Cjj  =  -  jf-S-  -  iVk^cosa  +  kySina)ds  (17) 

The  appropriate  boundary  condition  for  this  case  is  given  by 

2 

kx\|/2  +  kyXi/j  + - 2 — (k^cosa  +  kySina)  =  0  (18) 

yM^  P« 

Surface  Parameterization 

In  shape  design,  the  best  representation  of  the  surface  for  design  problems  remains  an 
open  issue.  In  the  current  study,  the  geometries  are  modeled  with  B  splines,  which  offer 
great  flexibility  in  the  definition  of  the  surfaces.  By  varying  the  polynomial  degree  and  the 
number  of  control  points,  a  wide  range  in  the  number  of  design  variables  and  surface  fi¬ 
delity  can  be  obtained.  On  one  hand,  the  design  variables  can  be  made  to  correspond  to  the 
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individual  grid  points  on  the  surface  by  choosing  a  linear  polynomial  and  an  appropriate 
number  of  control  points.  Conversely,  a  single  polynomial  curve  of  degree  n  (known  as  a 
Bezier  curve)  can  be  used  to  describe  the  geometry  by  choosing  the  number  of  control 
points  to  be  n  +  1 .  In  addition,  through  the  knot  sequence  associated  with  the  spline,  sharp 
breaks  in  the  surface  such  as  those  that  occur  in  cove  regions  and  blunt  trailing  edges  can 
still  be  represented  in  a  single  curve. 

In  a  B-spline  representation,  the  x-  and  y-coordinates  of  the  surfaces  are  written  in  a 
parametric  form  as'** 


n+1 

x(t)  =  I  X.N.  j^(t)  (19a) 

i  =  1 


n+  1 

y(t)  =  X  YjN.  ]^(t)  (19b) 

i  =  1 

where  (x,  y)  are  the  Cartesian  coordinates  of  the  surface,  N-  ,  is  the  B-spline  basis 
function  of  order  k,  (X-,  Yj^)  are  the  coordinates  of  the  B-spline  control  polygon,  and 
n+1  is  the  total  number  of  control  points.  Notice  that  the  surface  description  with  Eqs. 
(19a)  and  (19b)  is  still  continuous. 

In  Fig.  1,  a  point  on  the  old  surface  is  assumed  to  move  to  the  new  surface  while  re¬ 
maining  at  a  fixed  value  of  t .  Consequently,  variations  in  the  basis  functions  need  not  be 
considered.  In  addition,  generality  is  maintained  for  the  surface  geometry  as  variations  are 
not  restricted  to  being  strictly  normal  to  the  existing  surface. 


new  surface 


^  -  #-*’  t 


Figure  1.  Movement  of  point  on  surface. 
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For  small  variations,  an  incremental  length  on  the  new  surface  can  be  written  as 


where 


ds'  =  (1  +  A)ds 


C^x  +  CyY 

X  +y 


Cx=  IXi  -s 


(22a) 


C,=  X  Yi 


(22b) 


Here,  Xi  and  Yj  are  variations  in  the  position  of  the  B-spline  control  points,  and  x  and  y 
are  derivatives  with  respect  to  t .  Because  a  given  point  on  both  the  old  and  new  surfaces 
is  at  a  fixed  value  of  t ,  the  coordinates  on  the  new  surface  can  be  written  as 


=  IXrN,k=  I(Xi  +Xi)Nu 


=  I(Yi 


(23b) 


Therefore,  the  variation  of  a  point  on  the  surface  is  given  by 


=  IXiNi^k 


y=  lYiNi^k 


(24b) 


Since  the  components  of  the  surface  normal  can  be  expressed  as 


K  = 


I .2  .  .2 

Vx  +y 


(25a) 


ky  = 


-X 


(25b) 


the  variations  in  the  surface  normals  can  be  derived  as 


* '  r ,  L 


FTi  .  di 
vx  *  '• 


\\  +  y  i  =  1 


Hi  3. 
Vx  +y  1  = 


(26a) 


^  '  ~lT=2.^,'dl  -Tf=2  l  Yi; 

a/x  +  y  1  =  1  a/x  +  y  1  =  1 


‘Jt 


(26b) 


By  using  Eqs.  (13),  and  (20)-(26b),  variations  of  integrals  that  involve  gky  can  be  written 
as 


aQgkydsj  =  Jgkyds-  X  Xi 


Similarly,  variations  of  integrals  involving  gk^  are  given  by 


(27a) 


(27b) 


For  cost  functions  such  as  (3b)  that  only  involve  flow  quantities,  a  similar  procedure 
yields 
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(27c) 


sfjgdsl  =  jgds 

r 

YY-f  N  r.— +  _ — _ + 

^  Yr  ay  ; 


y^dt 


r 


+  IXil 

r 


N: 


dx  -2  .2Wt 
ox  ^  ui 


a/x^  +  dt 


The  terms  in  Eqs.  (27a),  (27b),  and  {21cl  involving  g  are  eliminated  by  the  boundary 
conditions  along  with  the  terms  involving  Q  in  Eq.  (12a)  as  discussed  earlier.  The  sensi¬ 
tivity  derivatives  are  obtained  by  combining  the  remaining  terms  in  Eq.  (27a),  (27b),  or 
(27c)  with  the  last  term  in  Eq.  (12a).  For  example,  using  Eqs.  (24a)-(26b)  to  compute  the 
variations  in  the  coordinates  and  metric  terms,  the  sensitivity  derivatives  of  the  drag  coef¬ 
ficient  with  respect  to  each  B-spIine  control  point  are  given  by 


yM, 


.agi  dg: 


dN 


i,  k 


dt 


+ 


|Ni,  k)(^i  +  u'J'2  +  +  H'F4)dt 


(28a) 


3c.i 


dYi  yM, 
2* 


dN 


i,k 


-J(-Q 

r 


dt 


+  u'P2  +  v'F3  +  H'F4)dt 


(28b) 


where  gj  =  (p/p^  -  l)cosa  and  g2  =  (p/p„  -  l)sina.  For  cost  functions  such  as  lift 
or  moment  coefficients,  a  similar  procedure  is  followed.  When  Mach  number  and  angle  of 
attack  are  considered  as  design  variables,  variations  from  surface  integrals  in  the  far  field 
also  contribute  to  the  sensitivity  derivatives. 

Navier-Stokes 

In  this  section,  the  adjoint  equations  with  the  associated  boundary  conditions  and  the  ex¬ 
pressions  for  the  sensitivity  derivatives  are  derived  for  viscous  flows.  Only  steady  incom¬ 
pressible  viscous  flows  are  considered  in  this  paper  to  make  the  analysis  more  transparent. 

The  governing  equations,  with  the  artificial  compressibility  parameter  p ,  are  given  by 
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(29) 


2 


+  p)+|.(uv)  =  i 


d  /^8u  8v 

ayUy  ax 


d .  \.a.2.  ,  u^a  /^au  av^ 

5^(UV)+5^V  =  + 


The  cost  function  is  augmented  as  in  Eq.  (2)  with  the  flow  equations  as  constraints 
through  the  Lagrange  multipliers  'P  with  Q  =  {p,  ^  v}  ^ .  The  variation  in  1^  is  split  into 
the  contributions  from  the  inviscid  and  the  viscous  terms  and  ,  respectively. 
These  can  be  derived  as 


61" 


“  ^  +  JQ’'(a’'1Cx  +  B^j'Pds 


(30) 


and 


—61 

11 


vise 

R 


-J5 


laxl  ax  )  ayUy  ax  j) 


dQ 


Q. 


'  3'r32 

a  f^'^2 

layt 

2 _ 

.  ay  y 

"all  j) 

dQ 


J'J'2[kx23^u  +  kyQyU  + 

-  +  k^j^^u  +  ^vjjds 


+ 


Sx  J 


ds 


+ 


a'F3 


ds 


(31) 


Combining  the  field  integrals  in  Eqs.  (30)  and  (31)  and  setting  the  integrands  to  zero 
yields  the  following  adjoint  system: 


_ 

ay 


(32a) 
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where 


0 

avo 

9 

ax' 

.  3x 

J  ay 

Uy  ax  J 

a'F3^ 

2-^) 

dx' 

k3y  dx  J 

ayl 

By  J 

(32b) 


For  purposes  of  illustration,  the  boundary  conditions  for  the  adjoint  system  are  derived 
with  the  assumption  that  the  cost  function  is  the  drag  coefficient; 


=  J 


c„k,,cosa  +  c„k„sina 

p  A  F  J 


2|i, 

Re 


2kj^^  +  k 

*dx 


cos  a 


Re 


(33) 


Note  that  the  expression  for  drag  retains  the  components  from  the  viscous  stress  tensor.  It 
will  be  seen  that  this  is  a  requirement  for  obtaining  boundary  conditions  for  the  adjoint 
equations.  Using  Eqs.  (27a)  and  (27b),  the  variation  in  the  drag  coefficient  is  given  by 


flow  ~  r/  •  ^Gi  •  9N:  .  \ 

=  &r  +  SXi|(YNu5^'-XNu5j"-5j'’''G2jdt 

302  3Ni,k. 


(34a) 


dt 


where 
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G,  =  2psina-4A|lsi„a-2i(|  +  |)cosa  (34d) 

Expressing  the  velocities  on  the  new  surface  in  a  Taylor  series  and  noting  that  the  veloci¬ 
ties  on  the  old  and  new  surface  are  both  zero,  the  variations  in  the  velocity  components 
can  be  written  as 


u 


V 


du  ~  9u  - 

dv  ~  dv  - 


(35a) 

(35b) 


In  order  to  derive  the  boundary  conditions,  Eqs.  (30),  (31),  (32a),  and  (32b)  are  com¬ 
bined,  and  terms  that  involve  the  variations  in  the  velocity  gradients  and  p  are  eliminated. 
This  requires  that  the  following  relationships  hold: 

k^'P2  +  ky'E3-i-2kj^cosa-i-2kySina  =  0  (36a) 

-4kjjCosa-2'P2kj^  +  4kySina  +  2kySina  =  0  (36b) 


2k„cosa-2k„sina-'Fok.. -'P^k..  =  0 


r36r'l 


This  system  is  overdetermined  and  is  satisfied  by  the  choice 


^^2  =  -2  COS  a 

(37a) 

^^3  =  -2  sin  a 

(37b) 

The  variation  in  the  drag  coefficient  can  be  obtained  from  Eq.  (34a)  by  using  these  equa¬ 
tions  in  conjunction  with  Eqs.  (24a)  and  (24b). 

Without  the  inclusion  of  the  full  stress  tensor  in  the  cost  function,  it  is  not  possible  to  ob¬ 
tain  a  consistent  set  of  boundary  conditions  for  4^2  and  4^3 .  Generally,  suitable  cost  func¬ 
tions  are  composed  of  terms  that  will  appropriately  balance  the  boundary  terms  from  the 
residuals.  In  particular,  cost  functions  such  as  lift,  drag,  and  pitching  moment  are  admissi¬ 
ble.  It  is  not  immediately  obvious  that  the  specification  of  a  pressure  distribution  is  allow¬ 
able  because  of  the  absence  of  viscous  terms  in  the  cost  function.  However,  a  suitable  cost 
function  can  be  obtained  by  first  replacing  the  pressure  term  in  the  stress  tensor  by  the  dif¬ 
ference  between  the  current  and  the  desired  pressure  coefficient  ACp .  This  is  then  premul¬ 
tiplied  by  the  surface  normal  scaled  by  this  difference  in  Cp  and  postmultiplied  by  the  sur¬ 
face  normal.  In  nondimensionalized  variables,  the  resulting  expression  is  given  by 


Ic  =  J{kxACp,kyACp} 

r 


2  Re  dx 
Revdy  dxj 


ji.  pu  dv^  ' 
Revdy  dxj 


2  Re  dy 


(38) 


After  expansion,  Eq.  (38)  can  be  rewritten  as 


-  r^ACpk,. 

Re  P 


2k„^  +  k 

*dx 


-f- 


^*^ydy  ^  ^ 


ds 


(39) 


This  equation  can  be  recast  in  terms  of  the  velocity  gradient  normal  to  the  boundary  as 


I 


C 


Re  dn  j 


ds 


(40) 


where  u„  is  the  normal  velocity  component  and  n  is  the  surface  normal  direction.  The  ve¬ 
locity  gradient  term  in  this  equation  is  zero  by  the  continuity  equation,  so  that  the  cost 
function  in  Eq.  (39)  corresponds  to  specifying  a  pressure  distribution.  However,  all  the 
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terms  in  Eq.  (39)  are  required  for  the  derivation  of  the  boundaiy  conditions  for  the  adjoint 
equations.  The  final  boundary  conditions  on  ^2  'P3  for  specifying  a  pressure  distri¬ 

bution  are  given  by 


'*'2  =  -2k,(Cp-c*) 

(41a) 

■Fa  =  -2k^(c^-c‘) 

(41b) 

The  continuous  adjoint  formulation  for  Navier-Stokes  equations  described  in  this  sec¬ 
tion  poses  a  problem  in  the  evaluation  of  the  sensitivity  derivatives.  The  evaluation  of 
these  derivatives  requires  second  derivatives  of  the  velocity  components  because  Gj  and 
G2  involve  velocity  gradients  that  are  further  differentiated  in  Eq.  (34a).  Recall  that  these 
terms  arise  from  expressing  the  cost  function  on  the  new  surface  in  a  Taylor  series  expan¬ 
sion  about  the  old  surface.  In  the  present  work,  because  the  flow  solver  is  only  second- 
order  accurate,  pointwise  second  derivatives  are  inconsistent  in  general.  An  accurate  eval¬ 
uation  of  second  derivatives  would  require  the  flow  solver  to  be  at  least  third-order  accu¬ 
rate.  If  a  mapping  is  employed,  as  is  possible  with  structured  grids,  the  surface  remains  at 
a  constant  coordinate  line,  and  this  problem  does  not  occur. 

In  a  discrete  adjoint  approach,  the  restriction  on  defining  a  suitable  cost  function  and  the 
need  for  second  derivatives  are  eliminated,  as  will  be  shown.  However,  the  full  implica¬ 
tion  of  designing  for  cost  functions  in  a  discrete  framework  for  which  boundary  condi¬ 
tions  are  not  obtainable  in  the  continuous  case  is  not  clear  at  this  time. 

Discretization 


Flow  equations 

The  discretization  of  the  flow  equations  is  first  addressed  since  it  has  implications  for 
the  discretization  of  the  adjoint  equations.  The  discretization  of  the  compressible  inviscid 
equations  is  given  first;  a  similar  procedure  is  used  to  discretize  the  inviscid  contributions 
for  the  incompressible  Navier-Stokes  equations.  The  equations  represent  a  system  of  con¬ 
servation  laws  for  a  control  volume  that  relates  the  rate  of  change  of  a  vector  of  state  vari¬ 
ables  Q  to  the  flux  through  the  volume  surface.  The  equations  are  written  in  integral  form 
as 


^jQdQ  +  ^F(Q,  n)ds  =  0  (42) 

where  for  compressible  flows  Q  =  [p,  pu,  pv,  E]^  and  F(Q,  n)  is  the  flux  of  mass,  mo¬ 
mentum,  and  energy  through  the  control  volume.  In  these  equations,  n  is  the  vector  nor¬ 
mal  to  the  boundary,  p  is  the  density,  u  and  v  are  the  Cartesian  velocity  components,  and 
E  is  the  total  energy  per  unit  volume.  These  equations  are  closed  by  the  equation  of  state 
for  a  perfect  gas. 
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In  discretizing  Eq.  (42),  the  variables  are  stored  at  the  vertices  of  a  triangular  mesh.  The 
control  volumes  are  defined  by  the  median  dual.  The  discrete  form  of  Eq.  (42)  for  vertex  i, 
with  an  associated  control  volume  Qj ,  is  given  by 

^jQd£2+  2;F,jl,j  =  0  ,43, 

j€Ni 

where  Fjj  is  the  numerical  flux  that  approximates  the  normal  flux  through  the  control-vol¬ 
ume  edge  dual  to  the  triangle  edge  that  joins  nodes  i  and  j ,  Ijj  is  the  length  of  the  dual 
edge,  and  Nj  is  the  set  of  vertex  neighbors  of  i .  The  numerical  fluxes  are  computed  by 
using  a  Roe-type  approximate  Riemann  solver:^'' 

Fy  =  |[F(Qj;n)-HF(Qj;n)-|A(Qr,Q,;n)|(Qr-Q,)]  (44) 

where  A  is  the  Jacobian  matrix  evaluated  at  the  Roe  state,  and  and  Q,  are  the  depen¬ 
dent  variables  on  the  right  and  left  boundaries  of  the  control  volume  face  which  are  ob¬ 
tained  by  extrapolation: 


Q,  +  fVQ. 

(fj-ri) 

(45a) 

Qj  +  fVQ. 

•(^i-rj) 

(45b) 

where  (p  =  0  for  first-order  discretization,  (p  =  1  for  second-order  discretization,  and  Fj 
and  Fj  are  the  position  vectors  of  nodes  i  and  j ,  respectively.  Note  that  the  definition  of 
the  fluxes  in  Eq.  (44)  is  different  from  a  standard  Riemann  solver  in  that  the  unsplit  fluxes 
are  evaluated  by  using  data  at  the  nodes  Qj  and  Qj  instead  of  data  at  the  extrapolated 
states  Qj.  and  Qj .  This  discretization  remains  second-order  accurate  and  has  the  benefit 
that  the  only  term  that  involves  data  other  than  at  the  immediate  neighbors  occurs  through 
the  dissipation.  This  enables  a  discretization  of  the  continuous  adjoint  equations  to  be  eas¬ 
ily  obtained  that  is  identical  to  the  discrete  adjoint  approach,  except  for  small  differences 
that  arise  from  the  higher  order  dissipation. 

For  computing  the  viscous  contributions  to  the  residual,  a  finite-volume  scheme  is  used 
that  is  equivalent  to  a  Galerkin  discretization  with  linear  basis  functions.  On  triangular 
grids,  this  discretization  only  requires  data  at  the  immediate  neighboring  nodes. 

Adjoint  equations 

The  adjoint  equations  can,  in  principle,  be  discretized  by  any  stable  and  consistent 
method.  However,  insufficient  grid  resolution  may  result  in  poor  accuracy  of  the  sensitiv¬ 
ity  derivatives  in  that  they  do  not  agree  with  those  obtained  by  finite  differences.  Inaccu¬ 
rate  sensitivity  derivatives  may  lead  to  failure  in  the  optimization  process.^’  Sensitivity  de- 
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rivatives  that  agree  with  finite-difference  gradients  can  be  obtained  regardless  of  grid  size 
by  making  the  equations  discretely  adjoint  to  the  discretized  flow  equations.  However, 
achieving  this  for  higher  order  discretizations  can  be  an  onerous  task.  In  the  present  work, 
the  discretization  is  been  derived  with  strong  guidance  from  a  discrete  adjoint  formulation 
so  that  for  first-order  accuracy,  a  direct  correspondence  with  a  discrete  adjoint  approach  is 
achieved.  Higher  order  accuracy  for  the  discretization  of  the  inviscid  terms  is  obtained 
through  the  use  of  extrapolation  of  the  costate  variables. 

The  discretization  of  the  adjoint  equation  is  performed  by  adding  a  time  derivative  to 
Eq.  (8)  and  using  a  finite-volume  type  of  method  similar  to  that  used  for  the  flow  solver. 
In  this  context,  Eq.  (8)  is  integrated  over  control  volumes,  where  the  matrices  are  taken 
outside  the  integrand  and  are  evaluated  using  nodal  point  values  of  the  dependent  vari¬ 
ables: 


(46) 


The  volume  integrals  are  converted  to  surface  integrals  over  each  of  the  control  volumes, 
and  the  values  of  the  costate  variables  on  the  boundaries  are  obtained  by  using  upwind 
type  formulas: 


(47) 


where  the  extrapolated  costate  variables  and  'Pj  are  obtained  by  using  formulas  that 
are  similar  to  Eqs.  (45a)  and  (45b).  The  data  used  for  evaluation  of  the  matrices  and  the 
formulas  used  for  obtaining  the  costate  variables  on  the  faces  of  the  control  volumes  have 
been  chosen  so  that  a  discrete  adjoint  formulation  is  obtained  for  first-order  spatial  accu¬ 
racy.  The  resulting  discretization  of  the  inviscid  contributions  may  be  written  as  follows: 


dt 


£2;  j  €  Nj 


(48) 


The  numerical  flux,  ,  used  in  calculating  the  residual  for  the  control  volume  that  sur¬ 
rounds  node  i ,  is  given  by 


:,j  =  5[A(Qi;n)Vi  +  'rj)  +  (||-)V-‘P,)' 


(49) 


where  O  =  |A(Qr,  Qj;n)|(Qj  -  Qj) .  Note  that  Gy  ^  -Gj; . 

On  solid  boundaries,  the  flux  along  the  wall  for  closing  off  the  surface  integral  around 
node  i  is  given  by 


Gjlj  =  k^Aj^'F  +  kyB/'^  =  aJ^'F 


(50) 
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For  enforcing  the  boundary  conditions  on  the  costate  variables,  a  weak  formulation  is  used 
in  which  the  fluxes  are  modified  appropriately  to  reflect  the  imposition  of  the  boundary 
conditions.  Numerically,  the  Jacobian  matrix  in  Eq.  (50)  is  evaluated  without  explicitly 
enforcing  the  boundary  condition  on  the  flow  variables  that  no  flow  is  allowed  through  the 
surface.  In  this  way,  the  contributions  from  the  fluxes  in  the  interior  in  conjunction  with 
the  boundary  flux  in  Eq.  (50)  combine  so  that  the  resulting  discretization  corresponds  with 
that  from  a  discrete  adjoint  approach. 

Note  that  in  Eq.  (49)  the  linearization  of  is  somewhat  cumbersome  but  has  been  pre¬ 
viously  derived  (see  for  example  Ref.  4).  A  simpler  equation  can  be  obtained  by  employ¬ 
ing  the  approximate  linearization  of  <1>  as 

i = lAi 

This  equation  is  less  complicated  than  the  full  linearization  and  only  differs  from  the  exact 
linearization  in  proportion  to  (Q^  -  Q]) ,  but  numerical  experiments  have  indicated  that  on 
very  coarse  grids  some  of  the  sensitivity  derivatives  are  of  poor  accuracy  compared  with 
finite-difference  derivatives.  Although  these  errors  decrease  as  the  grid  resolution  in¬ 
creases,  the  full  linearization  is  used  in  the  current  work. 

Since  the  viscous  equations  used  in  the  current  study  are  for  incompressible  flow,  the 
corresponding  terms  in  the  adjoint  equations  (Eq.  (32a))  have  the  same  form  and  are  there¬ 
fore  discretized  in  the  same  manner.  The  Dirichlet  boundary  conditions  for  Vj/j  and  Xj/j 
are  strongly  enforced  with  the  same  technique  used  to  set  the  velocities  to  zero  in  the  flow 
solver.  In  the  implicit  solver,  this  is  achieved  by  zeroing  the  off-diagonal  elements  in  the 
rows  of  the  matrix  that  correspond  to  boundary  nodes,  as  well  as  the  appropriate  terms  on 
the  right-hand  side. 

For  viscous  flows,  a  direct  correspondence  with  a  discrete  adjoint  formulation  is  not 
achieved  near  solid  boundaries.  This  is  easily  seen  by  examining  the  resulting  matrix 
stmctures  from  both  approaches  for  a  small  mesh  shown  in  Fig.  2,  where  it  is  assumed  that 
nodes  1 , 3,  and  5  lie  on  a  solid  wall. 


4  2 


Figure  2.  Sample  mesh. 
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In  the  discrete  adjoint  approach,  the  augmented  cost  function  is  given  by 


I(Q,  D,  'P,  X(D))  =  I,(Q,  D)  +  'P'^R(Q,  D,  X((D)))  (52) 

where  R  is  the  vector  of  discrete  residuals  and,  thus,  depends  explicitly  on  the  grid-point 
locations  X .  Taking  variations  of  Eq.  (52)  and  regrouping  terms  yields  the  adjoint  equa¬ 
tion 


r  IT 

aR 


.^Q 


ai 


c 


{'1'}  +  ^^!-  =  0 


The  variation  in  the  cost  function  is  then  given  by 


61  = 


c .  ^TraR  .  aRax>i 
aD  Ud  axaDj 


D 


(53) 


(54) 


In  these  equations,  it  is  understood  that  the  linearization  of  the  residual  includes  the  full 
effects  of  the  boundary  conditions.  Here,  aX/dD  represents  the  sensitivity  of  the  interior 
grid  points  to  changes  in  the  design  variables.  In  the  continuous  adjoint  formulation  de¬ 
scribed  earlier,  no  counterpart  to  this  term  exists.  The  determination  of  grid  sensitivities  is 
dependent  on  the  methodology  used  to  restructure  the  mesh.  Neglecting  these  terms  is 
equivalent  to  freezing  the  interior  grid  points,  regardless  of  changes  in  the  surface  geome¬ 
try.  Nevertheless,  in  a  second-order-accurate  scheme,  the  dR/dD  term  in  Eq.  (54)  ac¬ 
counts  for  changes  in  the  residuals  at  the  nodes  immediately  adjacent  to  the  surface,  as 
well  as  at  the  second  nearest  neighbors. 

A  diagram  of  the  matrix  structure  associated  with  the  configuration  of  nodes  in  Fig.  2  is 
shown  in  Fig.  3  for  the  discrete  adjoint  approach.  The  matrix  structure  for  the  continuous 
adjoint  approach  is  shown  in  Fig.  4.  In  these  figures,  the  solid  circles  represent  the  non¬ 
zero  entries  in  the  matrices.  Note  that  in  both  figures,  a  first-order  discretization  of  the  in- 
viscid  terms  is  assumed  so  that  the  stencil  only  involves  the  nearest  neighbors. 

Comparing  Figs.  3  and  4,  it  is  seen  that  the  matrix  structures  are  significantly  different. 
This  is  due  to  the  strong  enforcement  of  the  no-slip  condition  in  the  flow  solver,  which 
leads  to  zeros  in  the  columns  of  the  adjoint  system.  For  the  continuous  case,  explicit  en¬ 
forcement  of  the  boundary  condition  on  'Pj  and  ^^3  leads  to  zeros  along  rows.  Of  partic¬ 
ular  interest  in  the  discrete  adjoint  case  is  that  because  of  the  zeros  in  the  columns,  the  so¬ 
lution  of  the  costate  variables  in  the  interior  of  the  mesh  does  not  depend  on  the  values  of 
^2  and  ^3  at  the  boundary.  Furthermore,  because  the  residual  equation  for  the  flow 
solver  at  these  points  is  replaced  by  a  Dirichlet  condition  on  the  velocities,  the  residual 
does  not  depend  on  the  design  variables  so  that  8R/3D  =  0 .  Therefore,  there  is  no  con¬ 
tribution  to  the  sensitivity  derivatives  in  Eq.  (54)  from  these  terms.  The  result  is  that  in  the 
discrete  adjoint  case  the  values  of  'P2  and  ^^3  on  the  boundary  are  completely  arbitrary 
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and  have  no  effect  on  the  sensitivity  derivatives.  This  has  been  verified  by  numerical  ex¬ 
periments. 


Columns  associated  with  Dirichlet  rows 


Figure  3.  Matrix  structure  for  discrete  adjoint  approach. 


Row  with  Dirichlet 
boundary  conditions 


Row  with  Dirichlet 
boundary  conditions 


Row  with  Dirichlet 
boundary  conditions 

Figure  4.  Matrix  structure  for  continuous  adjoint  approach. 

In  light  of  the  discussion  above,  it  is  of  interest  to  compare  the  values  of  the  costate  vari¬ 
ables  that  are  obtained  from  both  the  continuous  and  the  discrete  adjoint  formulations  for 
a  viscous  flow.  In  Fig.  5,  profiles  of  '1^2  as  a  function  of  the  distance  from  the  body  are 
shown  for  a  case  in  which  the  cost  function  is  the  drag  of  an  airfoil  and  the  location  of  the 
profile  is  taken  to  be  at  the  midchord  of  the  airfoil  on  the  upper  surface.  In  the  figure,  the 
values  of  'Fj  agree  well  away  from  the  body.  Near  the  boundary,  however,  the  costate 
variables  from  the  continuous  and  discrete  formulations  do  not  agree.  As  the  mesh  is  re¬ 
fined,  the  distance  from  the  surface  of  the  airfoil  in  which  these  discrepancies  occur  de- 
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creases.  Thus,  one  would  expect  that  in  the  limit  of  vanishing  mesh  size,  the  two  ap¬ 
proaches  would  agree. 


Coarse 


-  Discrete 
Continuous 

Medium  Fine 


Xg 


Xg 


Figure  5.  Profiles  of  obtained  from  both  discrete  and  continuous  adjoint  formulations. 

In  Eq.  (53),  'F  can  be  determined  provided  that  dR/9Q  is  nonsingular  irrespective  of 
the  cost  function.  Also,  no  difficulty  is  encountered  in  determining  the  sensitivity  deriva¬ 
tives  with  Eq.  (54).  In  particular,  note  that  this  equation  does  not  require  explicit  calcula¬ 
tion  of  second  derivatives.  Therefore,  for  viscous  flows,  a  discrete  approach  is  used  in  the 
current  study,  except  that  higher  order  accuracy  for  the  inviscid  terms  is  achieved  by  using 
the  continuous  approach  described  in  the  inviscid  section.  The  implementation  of  this  ap¬ 
proach  does  not  entail  much  additional  effort  because  the  inviscid  terms  are  already  dis¬ 
cretely  adjoint  for  first-order  accuracy,  and  the  viscous  terms  only  involve  the  nearest 
neighbors.  The  accuracy  of  the  derivatives  using  this  approach  is  comparable  to  that  ob¬ 
tained  for  inviscid  flows.  For  first-order  accuracy,  the  resulting  method  is  identical  to  the 
standard  discrete  adjoint  approach. 

Solution  procedures 

For  the  flow  equations,  an  implicit  solution  method  with  multigrid  acceleration  is  used. 
Details  may  be  found  in  Refs.  1,  2,  and  9.  The  discretized  equations  for  the  costate  vari¬ 
ables  in  the  absence  of  the  time  derivative  represent  a  linear  system  that  can  be  solved  by 
using  a  technique  such  as  preconditioned  GMRES.^^  Alternatively,  by  retaining  the  time 
derivative,  the  equations  can  be  solved  to  steady  state  by  using  a  time-marching  proce¬ 
dure.  In  the  present  work,  the  time  term  is  included  and  a  multigrid  procedure  is  used  with 
preconditioned  GMRES  as  a  smoother.  The  preconditioning  is  accomplished  using  an  in¬ 
complete  lower/upper  (LU)  decomposition  with  no  fill-in.  The  motivation  for  retaining 
the  time  term  is  that  this  approach  often  converges  in  situations  for  which  the  GMRES 
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procedure  might  otherwise  “stall.”  Note  that  because  the  equations  are  linear,  the  matrix- 
vector  products  are  easily  formed  by  simply  passing  the  vector  to  the  residual  routine  in 
place  of  the  costate  variables.  By  forming  the  matrix-vector  products  in  this  way,  the  larg¬ 
est  contribution  to  memory  requirements  is  through  the  preconditioner  so  that  the  result¬ 
ing  scheme  requires  roughly  the  same  amount  of  memory  as  the  flow  solver. 

Grid  generation  and  mesh  movement 

The  unstructured  meshes  used  in  this  work  are  generated  using  the  software  package  de¬ 
scribed  in  Ref.  28.  This  employs  an  advancing  front  type  of  method  that  generates  good 
quality  grids  for  both  inviscid  and  viscous  calculations. 

For  shape  optimization,  the  design  is  carried  out  in  a  domain  that  changes  during  the  de¬ 
sign  cycle  as  the  shape  of  the  boundary  changes.  Therefore,  the  existing  grid  is  modified 
in  order  to  conform  to  the  changing  domain. 

For  inviscid  flows,  the  strategy  outlined  in  Ref.  42  is  used  to  restructure  the  mesh  in  re¬ 
sponse  to  the  changes  in  the  surface  shape.  The  tension-spring  analogy  is  employed  to 
allow  the  field  grid  points  to  respond  to  the  displacements  of  the  points  on  the  surface.  The 
following  linear  system  of  equations  is  solved  with  a  Jacobi  iteration  strategy: 

j€  Nj 

where  AXj  and  Axj  are  the  displacements  from  the  initial  position  for  nodes  i  and  j.  The 
spring  stiffness  is  assumed  to  be  Ijj”  ,  where  Ijj  is  the  length  of  the  edge  that  joins 
nodes  i  and  j .  Note  that  by  using  Eq.  (55),  the  mesh  remains  unchanged  when  the  surface 
is  held  fixed.  When  the  boundary  shape  changes  during  the  design  cycle,  this  method  does 
not  guarantee  that  the  grid  lines  will  not  cross.  An  improvement  is  to  make  the  spring  sys¬ 
tem  nonlinear  (i.e.,  the  shape  change  is  decomposed  into  smaller  steps,  and  the  procedure 
is  repeated  at  each  step).  Also,  in  order  to  maintain  good  mesh  quality  throughout  the  de¬ 
sign  cycle,  the  edges  are  reconnected  according  to  either  a  Delaunay  criterion  or  by  locally 
minimizing  maximum  angles  (min-max). 


Figure  6.  Methodology  for  mesh  movement  for  viscous  grids. 
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For  grids  with  high-aspect-ratio  cells,  the  inviscid  strategy  fails  for  a  number  of  reasons. 
The  spring  analogy  typically  results  in  invalid  grids  with  crossing  of  grid  lines.  In  addi¬ 
tion,  both  the  Delaunay  and  min-max  criteria  often  result  in  nodes  with  large  connectivi¬ 
ties.  Therefore,  the  grid-movement  scheme  is  modified  to  deal  with  Navier-Stokes  grids. 
The  Delaunay  criterion  is  replaced  by  the  min-max  criterion  where  the  swapping  is  only 
carried  out  if  the  maximum  angle  exceeds  a  specified  angle  (set  to  150° ).  The  distance  to 
the  wall  for  each  node  in  the  mesh  is  first  computed.  When  the  points  on  the  surface  are 
displaced,  the  field  points  move  in  response,  as  shown  in  Fig.  6(a).  Here,  AB  is  an  edge  on 
the  surface  of  the  body.  Nodes  A  and  B  move  to  A’  and  B’,  respectively.  For  the  field 
point  X,  the  nearest  point  on  edge  AB  is  denoted  by  C.  Given  vectors  AA’  and  BB’,  the 
vector  CC’  is  obtained  by  linear  interpolation.  The  field  point  X  moves  to  X’  such  that 
XX’  is  equal  and  parallel  to  CC’.  In  order  to  contain  the  effect  of  grid  movement  to  a  spec¬ 
ified  region,  XX’  is  multiplied  by  an  exponential  factor  that  decays  from  unity  at  the  sur¬ 
face  to  nearly  zero  at  a  specified  cut-off  distance.  This  technique,  in  combination  with 
edge  swapping,  allows  for  large  changes  in  body  shapes  even  when  highly  stretched  grids 
are  used.  However,  the  grids  tend  to  lose  orthogonality  near  the  surface  when  large 
changes  occur  in  the  surface  geometry.  To  improve  orthogonality  near  the  surface,  the 
method  described  above  is  replaced  by  the  one  shown  in  Fig.  6(b)  within  a  specified  dis¬ 
tance  to  the  wall.  In  this  technique,  CC’  is  obtained  as  before,  but  C’X’  remains  orthogo¬ 
nal  to  A’B’  and  the  normal  distance  d  is  maintained.  It  is  also  desirable  to  revert  to  the  in¬ 
viscid  algorithm  in  regions  where  the  grid  is  not  highly  stretched.  Therefore,  outside 
another  specified  distance  from  the  wall,  the  inviscid  algorithm  is  employed.  Thus,  the 
final  scheme  is  a  blending  of  all  three  methods.  This  scheme  has  been  found  to  be  effec¬ 
tive  in  dealing  with  Navier-Stokes  grids,  even  for  large-scale  changes  in  surface  shape, 
and  is  reasonably  insensitive  to  the  cutoff  distances  provided  that  the  region  in  which  or¬ 
thogonality  is  maintained  is  restricted  to  the  immediate  vicinity  of  the  wall.  Unless  the  dis¬ 
placements  of  the  surfaces  are  large,  the  last  step  can  be  skipped. 

The  technique  described  above  is  demonstrated  in  Fig.  7  for  a  Navier-Stokes  2rid  about 
an  airfoil.  The  grid  contains  26949  nodes,  and  the  spacing  at  the  wall  is  2  x  10^  relative 
to  the  chord.  In  this  figure,  the  nose  of  the  geometry  is  distorted  by  moving  one  of  the 
Bezier  control  points  in  this  region.  Although  the  geometry  is  significantly  altered,  a  valid 
mesh  results,  which  maintains  good  quality  as  well  as  orthogonality  near  the  surface.  It 
should  be  pointed  out  that  for  multielement  configurations,  the  procedure  described  may 
fail  for  large  relative  displacements  of  the  elements  because  the  cutoff  regions  that  may  be 
initially  distinct  could  “collide.”  Further  work  is  necessary  in  this  area. 
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Figure  7.  Example  of  mesh  movement  for  viscous  mesh. 

Optimizer 

The  optimizer  used  in  the  current  study  is  KSOPT/®  which  uses  a  quasi-Newton  method 
to  determine  the  search  directions  and  a  polynomial  line  search  technique  to  determine  the 
step  length  in  the  descent  direction.  This  code  has  been  chosen  because  it  is  capable  of 
multipoint  design  and  can  handle  both  equality  and  inequality  constraints.  In  addition, 
upper  and  lower  bounds  can  be  placed  on  design  variables;  this  is  the  method  that  is  cur¬ 
rently  used  to  enforce  the  geometric  constraints  necessary  to  maintain  a  viable  geometry 
throughout  the  design  cycle. 


Results 


Accuracy  of  derivatives 

To  assess  the  accuracy  of  derivatives,  an  isolated  transonic  airfoil  and  a  subsonic  mul¬ 
tielement  airfoil  (where  interaction  between  the  elements  occurs  through  the  flow  field) 
are  studied.  For  the  first  test,  a  single  12th-order  Bezier  curve  is  used  to  approximate  an 
NACA  0012  airfoil,  with  only  13  control  points.  In  the  experiment  that  follows,  a  grid 
with  4770  nodes  is  generated,  with  128  grid  points  on  the  surface  of  the  airfoil.  The  cost 
function  is  the  lift  coefficient,  and  derivatives  with  respect  to  the  Bezier  control  points  are 
obtained  using  the  continuous  adjoint  method  and  are  compared  with  those  from  finite  dif¬ 
ferences.  The  Mach  number  for  this  case  is  0.75,  and  the  angle  of  attack  is  1.25° .  The  re¬ 
sulting  pressure  distribution  is  shown  in  Fig.  8  and  exhibits  a  shock  on  the  upper  surface 
of  the  airfoil. 
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Figure  8.  Pressure  distribution  for  NAC A  0012  with  =  0.75  and  a  =  1.25°. 

The  sensitivity  derivatives  of  the  lift  with  respect  to  the  y  position  of  the  individual  con¬ 
trol  points  are  shown  in  Figs.  9  and  10  using  the  continuous  adjoint  approach.  In  these  fig¬ 
ures,  the  derivatives  are  obtained  by  using  the  second-order  formulation  for  both  the  flow 
solver  and  the  adjoint  equations.  The  corresponding  derivatives  for  first-order  accuracy 
are  not  shown  because  the  first-order  scheme  has  been  verified  to  be  discretely  adjoint  to 
the  flow  solver  in  this  case.  In  Fig.  9,  the  derivatives  at  the  first  and  last  control  points 
(numbers  1  and  13)  correspond  to  those  at  the  trailing  edge.  Although  the  derivatives  of 
the  control  points  at  the  trailing  edge  are  available  from  the  adjoint  approach,  the  corre¬ 
sponding  finite-difference  derivatives  are  not  obtained  because  the  geometry  would  “sep¬ 
arate”  at  the  trailing  edge.  Instead,  the  grid  point  at  the  trailing  edge  is  perturbed,  and  the 
resulting  derivative  is  compared  with  the  sum  of  the  derivatives  at  this  location  from  the 
adjoint  approach.  A  close-up  view  of  the  derivatives  away  from  the  trailing  edge  is  shown 
in  Fig.  10.  The  figures  indicate  that  the  derivatives  are  fairly  accurate;  the  largest  discrep¬ 
ancy  between  the  adjoint  and  the  finite-difference  derivatives  is  less  than  5  percent.  Note 
that  in  this  study  specification  of  the  costate  variables  as  a  boundary  condition  across  dis¬ 
continuities  in  the  field,  as  suggested  in  Refs.  17  and  18,  is  not  done  with  no  apparent  deg¬ 
radation  in  accuracy. 
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Figure  9.  Comparison  of  derivatives  obtained  using  adjoint  approach  with  finite 
differences  for  NAC A  0012  with  =  0.75  and  a  =  1.25°. 


Control  Point 

Figure  10.  Close-up  view  of  Fig.  9 

As  mentioned  previously,  unless  the  adjoint  equations  are  discretized  appropriately,  the 
resulting  derivatives  may  exhibit  inaccuracies  when  compared  with  finite-difference  gra¬ 
dients  on  coarse  grids.  To  study  this  aspect  more  closely,  a  two-element  airfoil  is  consid¬ 
ered  for  which  the  surface  of  each  element  is  represented  with  a  third-order  B-spline  with 
31  control  points.  The  cost  function  is  the  lift  coefficient,  and  the  derivatives  with  respect 
to  the  design  variables  on  the  aft  element  are  computed  with  both  methods  on  a  set  of  four 
sequentially  finer  grids.  These  grids,  denoted  as  grids  4,  3, 2,  and  1,  consist  of  1 103, 3030, 
9591,  and  18,392  nodes,  respectively;  of  these  nodes,  88,  176, 352,  and  704  lie  on  the  air- 
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foil  surfaces.  Obtaining  the  sensitivity  derivatives  with  central  finite-difference  formulas 
requires  58  flow-field  computations  for  each  grid.  For  the  adjoint  approach,  all  derivatives 
are  obtained  in  one  solution  of  the  adjoint  equations,  which  requires  roughly  the  same 
amount  of  work  as  one  solution  of  the  flow  field. 


CD 


CO 


Design  Variable  Design  Variable 


Figure  11.  Accuracy  of  derivatives  on  aft  element  of  a  two-element  airfoil;  third-order  B- 

spline  with  31  control  points 

In  Fig.  11,  finite-difference  derivatives  are  compared  with  those  obtained  using  the  for¬ 
mulas  for  the  continuous  adjoint  approach.  In  this  figure,  the  derivatives  in  the  immediate 
vicinity  of  the  trailing  edge  are  not  shown  so  that  the  derivatives  over  the  bulk  of  the  air¬ 
foil  can  be  examined  more  closely.  The  importance  of  the  derivatives  near  the  trailing 
edge  is  discussed  later  in  this  section.  In  addition,  derivatives  are  also  shown  from  a  “hy¬ 
brid”  approach  in  which  the  costate  variables  are  obtained  from  the  continuous  adjoint  ap¬ 
proach  and  are  subsequently  used  in  a  discrete  adjoint  framework  to  compute  the  sensitiv¬ 
ity  derivatives  by  using  Eq.  (54).  In  this  approach,  no  approximations  are  used  in  Eq.  (54), 
so  that  the  only  difference  between  the  hybrid  approach  and  a  purely  discrete  adjoint  ap¬ 
proach  stems  from  small  differences  in  obtaining  the  costate  variables  for  the  second- 
order  discretization.  Recall  that  in  the  continuous  adjoint  case,  no  sensitivities  that  result 
from  mesh  movement  appear  in  the  equations.  Therefore,  for  the  hybrid  approach,  the 
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mesh  sensitivities  have  been  neglected  in  order  to  examine  the  effects.  As  discussed  ear¬ 
lier,  when  the  surface  of  the  airfoil  is  perturbed,  the  residuals  are  affected  at  the  nodes  on 
the  body  as  well  as  at  their  first  and  second  nearest  neighboring  nodes  (the  second  nearest 
nodes  are  affected  through  the  gradient  computation).  Therefore,  the  residuals  at  these 
nodes  contribute  to  the  sensitivity-derivative  calculation  in  Eq.  (54).  However,  inclusion 
of  these  contributions  does  not  account  for  any  change  in  the  interior  residuals  caused  by 
the  possible  movement  of  interior  mesh  points.  The  situation  is  thus  equivalent  to  the  case 
for  which  the  surface  of  the  airfoil  is  modified  but  the  interior  of  the  mesh  is  held  fixed. 

As  seen  in  the  figure,  the  finite-difference  derivatives  are  nonsmooth  on  the  coarser 
grids  and  have  several  derivatives  of  negative  sign.  The  derivatives  obtained  with  the  hy¬ 
brid  approach  follow  an  almost  identical  pattern.  The  derivatives  obtained  from  the  con¬ 
tinuous  adjoint  approach  are  smoother  on  all  of  the  grids  and  remain  positive  over  the  en¬ 
tire  interval  shown.  Although  discrepancies  result  over  parts  of  the  airfoil,  the  derivatives 
calculated  with  all  three  methods  agree  as  the  grids  are  refined.  A  case  could  be  made  that 
the  continuous  adjoint  derivatives  are  “better”  because  the  signs  of  derivatives  are  always 
in  “correct”  agreement  with  those  from  the  finest  grid.  However,  when  designing  on  the 
coarser  grids,  this  could  cause  the  optimizer  to  fail  because  the  derivatives  do  not  accu¬ 
rately  represent  the  discrete  derivatives.^’  Conversely,  the  hybrid  approach  may  be  consid¬ 
ered  to  be  “better”  in  that  the  derivatives  agree  more  closely  with  those  obtained  from  fi¬ 
nite  differences  on  all  the  meshes.  Although  this  may  lead  to  successful  numerical 
optimization  on  all  grids,  the  resulting  geometry  may  be  quite  different  from  that  obtained 
with  a  finer  grid.  In  either  case,  a  suitably  refined  grid  must  be  employed  in  which  case 
neither  the  continuous  nor  the  discrete  approach  offers  a  significant  advantage  over  the 
other. 

In  Fig.  11,  the  discrepancies  in  the  derivatives  on  the  coarse  grids  stem  from  three 
sources.  These  include  the  fact  that  the  second-order  scheme  is  not  exactly  discretely  ad¬ 
joint  to  the  flow  equations  on  all  grids.  Also,  small  errors  in  the  finite-difference  calcula¬ 
tions  may  be  present  as  a  result  of  the  choice  of  step  size  which  was  not  optimized  for  ac¬ 
curacy  for  each  of  the  29  design  variables  although  a  reasonable  effort  was  made  to 
determine  acceptable  values.  In  addition,  the  derivatives  obtained  from  finite  differences 
include  the  effect  of  grid  sensitivities  because  the  interior  mesh  points  are  relaxed  each 
time  a  design  variable  is  perturbed  using  the  techniques  described  earlier.  As  mentioned 
previously,  neither  the  continuous  adjoint  nor  the  hybrid  approach  has  included  these  ef¬ 
fects  because  the  continuous  formulation  assumes  no  dependence  on  a  grid  and  the  hybrid 
formulation  has  neglected  these  contributions  for  this  test.  The  figure  shows  clearly  that  as 
the  grids  are  refined  the  derivatives  over  the  bulk  of  the  airfoil  approach  the  same  value 
regardless  of  the  methodology  used  to  obtain  them. 

Although  it  is  tempting  to  conclude  from  the  above  example  that  grid  sensitivities  do  not 
play  a  major  role  as  the  grid  is  refined,  this  conclusion  is  not  always  valid.  To  demon¬ 
strate,  a  simple  example  is  given  in  which  the  geometry  and  flow  conditions  are  held  fixed 
while  the  grid  is  allowed  to  change.  More  specifically,  the  relationship  between  the  airfoil 
surface  and  the  grid  is  changed.  The  role  of  the  grid  sensitivities  is  studied  by  considering 
the  derivative  of  the  lift  of  a  single  airfoil  due  to  a  vertical  translation. 

For  this  case,  an  NACA  0012  airfoil  at  a  free-stream  Mach  number  of  0.5  and  an  angle 
of  attack  of  2°  is  considered.  A  sequence  of  structured  C-type  grids  is  utilized  in  which 
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each  grid  represents  a  uniform  refinement  in  each  direction  over  the  previous  level.  Two 
structured-grid  codes^’-'*’  are  used,  in  addition  to  the  unstructured-grid  flow  solver.  For  the 
unstructured  flow  solver,  the  cells  in  the  structured  mesh  are  simply  divided  into  triangles. 
The  derivative  of  the  lift  with  respect  to  translation  of  the  airfoil  surface  in  the  y  direction 
is  computed  with  central  differences.  The  airfoil  surface  is  perturbed  a  small  amount,  and 
three  different  techniques  are  considered  for  modifying  the  interior  mesh: 

1 .  The  airfoil  surface  and  the  entire  mesh  are  shifted. 

2.  The  airfoil  surface  is  perturbed,  and  the  rest  of  the  mesh  remains  fixed. 

3.  The  airfoil  surface,  as  well  as  the  mesh  line  that  extends  from  the  trailing  edge  of  the 
airfoil  to  the  downstream  outer  boundary  are  perturbed,  and  the  rest  of  the  grid 
remains  fixed. 

For  the  case  in  which  the  airfoil  and  the  mesh  are  simultaneously  perturbed,  the  lift  does 
not  change  and  the  derivative  is  zero,  independent  of  the  mesh  size.  This  case  corresponds 
to  simply  a  shifting  of  the  origin  of  the  coordinate  system;  therefore,  no  calculations  are 
required.  The  importance  of  the  second  method  for  computing  the  finite-difference  deriv¬ 
ative  is  that  this  situation  corresponds  to  the  case  in  which  grid  sensitivities  are  ignored  in 
a  discrete  formulation.  This  correspondence  has  been  verified  using  the  derivatives  ob¬ 
tained  from  the  first-order  adjoint  code,  where  the  derivatives  are  obtained  by  using  the 
hybrid  methodology  and  the  grid  sensitivities  are  neglected.  The  third  method  is  chosen 
simply  for  demonstration  purposes.  Note  that  in  the  numerical  experiments  that  follow,  all 
results  are  converged  to  machine  zero  and  the  step  size  for  computing  the  finite-difference 
derivatives  has  been  varied  over  a  large  range  of  values  with  no  significant  changes  in  the 
results.  In  all  cases,  the  step  size  that  is  used  is  much  smaller  than  the  distance  from  the 
surface  of  the  airfoil  to  the  first  grid  line,  so  no  crossing  of  grid  lines  occurs. 

In  an  ideal  situation,  the  lift  of  a  single  airfoil  in  an  unbounded  flow  would  be  insensi¬ 
tive  to  a  vertical  change  in  the  coordinates  so  that  the  derivative  would  be  zero.  Numeri¬ 
cally,  however,  changes  may  occur  because  of  the  changing  location  of  the  airfoil  relative 
to  the  outer  boundary  and  because  of  possible  changes  in  the  grid.  In  the  case  where  the 
entire  grid  is  shifted,  the  derivative  of  lift  due  to  a  shift  in  the  y  location  of  the  surface  is 
zero.  By  shifting  all  grid  lines  except  the  one  at  the  outer  boundary,  the  derivatives  have 
been  found  to  remain  very  small  (0(  10"^) )  which  indicates  that  the  derivative  of  lift  due 
to  the  location  of  the  outer  boundary  is  small.  In  this  case,  the  changes  are  not  only  attrib¬ 
utable  to  the  changing  location  of  the  outer  boundary  but  also  to  some  small  grid  effects  at 
the  outer  boundary. 

In  Fig.  12,  the  sensitivity  derivatives  of  the  lift  with  respect  to  translation  of  the  airfoil  in 
the  y  direction  are  shown  for  methods  2  and  3  described  above.  As  seen  in  the  figure,  the 
derivatives  due  to  the  translation  of  the  airfoil  surface  depend  greatly  on  the  methodology 
used  to  modify  the  grid.  More  importantly,  these  derivatives  do  not  tend  to  zero  as  the 
mesh  is  refined  but  actually  increase  in  magnitude! 
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Figure  12.  Finite-difference  derivatives  of  lift  with  respect  to  vertical  shift  in  airfoil 
position  obtained  with  fixed  computational  grid. 

Computing  the  derivative  of  lift  with  respect  to  a  vertical  translation  corresponds  to  a 
simple  summation  of  the  derivatives  of  lift  with  respect  to  the  y  position  of  each  of  the  de¬ 
sign  variables.  For  example,  3cj/8Y  for  the  aft  element  of  the  airfoil  shown  in  Fig.  1 1 
can  be  computed  by  a  summation  of  the  individual  derivatives.  Although  the  individual 
derivatives  shown  in  Fig.  1 1  converge  as  the  mesh  is  refined,  the  derivatives  at  the  trailing 
edge  do  not.  In  Fig.  13,  the  individual  sensitivity  derivatives  that  are  obtained  with  the  hy¬ 
brid  approach  are  now  plotted  at  a  scale  so  that  the  derivatives  at  the  trailing  edge  can  be 
seen.  Whereas  the  derivatives  away  from  the  trailing  edge  converge  as  the  grid  is  refined 
(see  Fig.  1 1),  those  at  the  trailing  edge  of  the  airfoil  do  not  exhibit  the  same  level  of  con¬ 
vergence  and,  in  fact,  continually  change  as  the  grid  is  refined.  This  behavior  appears  to 
be  caused  by  the  singularity  at  the  trailing  edge  and  is  the  source  of  the  sensitivity  of  the 
derivative  to  the  manner  in  which  the  grid  is  treated. 


30 


Figure  13.  Sensitivity  derivatives. 


In  Fig.  14,  finite-difference  derivatives  similar  to  those  shown  in  Fig.  12  are  shown  for  a 
symmetric  Joukowski  airfoil  at  the  same  Mach  number  and  angle  of  attack  as  before.  For 
this  airfoil,  the  slope  of  the  upper  and  lower  surfaces  in  the  analytical  definition  are  identi¬ 
cal  at  the  trailing  edge,  and  the  effect  of  the  singularity  should  be  reduced.  These  deriva¬ 
tives  have  been  obtained  by  shifting  only  the  surface  of  the  grid,  as  in  method  2.  As  is 
clearly  seen  in  the  figure,  the  derivatives  are  much  smaller  in  magnitude  than  those  for  the 
NACA  0012  and  do  not  increase  in  magnitude  as  the  grid  is  refined. 

From  the  foregoing  discussion,  it  is  apparent  that  the  grid  sensitivities  near  the  trailing 
edge  of  the  airfoil  can  play  a  major  role  in  the  computation  of  the  derivatives  necessary  to 
position  airfoils  relative  to  one  another.  It  should  be  emphasized  that  during  an  actual  de¬ 
sign  the  grid  is  generally  “relaxed,”  so  that  the  original  relationships  between  the  grid 
points  are  more  or  less  intact,  and  that  the  effect  of  the  grid  would  be  much  less  pro¬ 
nounced  than  that  shown  above.  The  important  point  is  that  without  inclusion  of  the  grid 
sensitivities,  the  derivatives  obtained  would  correspond  to  the  case  above  in  which  the  in¬ 
terior  grid  is  held  fixed.  Because  the  derivatives  clearly  depend  on  the  manner  in  which 
the  mesh  and  the  geometry  interact,  this  factor  must  be  accounted  for  in  the  computations 
when  derivatives  are  needed  in  the  immediate  vicinity  of  the  trailing  edge.  Furthermore, 
the  errors  caused  by  failure  to  properly  account  for  these  terms  do  not  vanish  as  the  mesh 
is  refined.  However,  from  the  results  in  Fig.  11,  it  appears  that  grid  sensitivities  can  be 
safely  neglected  in  regions  away  from  the  trailing  edge,  provided  that  the  grid  is  suffi¬ 
ciently  refined. 
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Figure  14.  Translation  for  Joukowski  airfoil. 

Inviscid  Design  Examples 

An  example  of  shape  optimization  is  shown  below  in  which  drag  minimization  has  been 
performed  for  a  single  airfoil.  The  initial  geometry  is  an  NACA  0012  airfoil,  described  by 
a  third-order  B-spline  with  50  control  points.  The  grid  consists  of  4763  nodes  with  128 
nodes  on  the  airfoil  surface.  The  Mach  number  for  this  test  is  0.75,  and  the  initial  angle  of 
attack  is  2° .  For  this  case,  the  computed  lift  coefficient  is  0.4229,  with  a  corresponding 
drag  coefficient  of  0.0123.  For  this  design,  the  cost  associated  with  maintaining  the  cur¬ 
rent  lift  coefficient  is  combined  with  that  for  minimizing  the  drag: 

Ic(Q,  D)  =  ^(ci  -  c,*)^  +  10  X  l(Cd  -  (56) 

where  Cj*  is  the  desired  lift  coefficient  and  c^*  is  zero.  The  factor  of  10  associated  with 
the  cost  function  for  drag  is  chosen  so  that  the  contribution  from  each  cost  function  is  of 
the  same  order  of  magnitude.  The  design  variables  are  the  y-coordinates  of  the  control 
points  that  describe  the  airfoil,  except  those  at  the  trailing  edge,  which  remain  fixed.  The 
angle  of  attack  is  an  additional  design  variable  and  is  allowed  to  vary  in  order  to  maintain 
the  lift.  The  total  number  of  design  variables  for  this  case  is  49.  For  this  case,  the  continu¬ 
ous  adjoint  approach  is  used  instead  of  the  hybrid  approach. 

After  10  design  iterations,  the  lift  coefficient  is  0.4225,  which  is  in  close  agreement 
with  the  specified  lift  coefficient  of  0.4229.  The  drag  has  been  reduced  from  0.0123  to 
0.0016,  and  the  final  angle  of  attack  is  1.747° .  The  objective  function  and  the  root  mean 
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square  (rms)  of  the  sensitivity  derivatives  have  each  been  reduced  between  1  and  2  orders 
of  magnitude.  Note  that  these  gradients  are  not  the  projected  gradients  and  that  several 
side  constraints  are  active.  The  initial  and  final  pressure  distributions  are  shown  in  Fig.  15; 
the  corresponding  geometries  are  shown  in  Fig.  16. 
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Figure  15.  Initial  and  final  pressure  distribution  for  NACA  0012  design  case. 
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Figure  16.  Initial  and  final  airfoil  for  NACA  0012  design  case. 


The  next  case  is  that  of  a  two-element  airfoil  configuration  that  consists  of  two  airfoils 
in  which  the  top  airfoil  is  displaced  from  the  other  in  the  positive  y  direction  by  0.5  chord 
and  in  the  negative  x  direction  by  0.5  chord.  The  free-stream  Mach  number  is  0.60,  and 
the  angle  of  attack  is  0° .  A  sequence  of  three  grids  for  use  with  multigrid  acceleration  has 
been  generated  for  this  case.  The  finest  grid  consists  of  7974  points  and  is  shown  in  Fig. 
17. 


Figure  17.  Initial  configuration  for  two-element  test  case. 

For  this  case,  the  objective  is  to  modify  the  shape  of  the  aft  airfoil  in  order  to  achieve  a 
desired  pressure  distribution  on  the  front  airfoil.  The  desired  pressure  distribution  has 
been  obtained  from  analysis  of  the  initial  configuration,  with  the  shape  of  the  aft  airfoil 
modified.  Although  this  test  case  is  somewhat  fabricated,  it  demonstrates  flexibility  that 
may  be  difficult  to  achieve  with  inverse  methods  in  which  the  interaction  between  ele¬ 
ments  is  not  taken  into  account. 

Pressure  contours  for  the  initial  flow  field  are  shown  in  Fig.  18(a);  the  corresponding 
contours  of  ^2  shown  in  Fig.  18(b).  The  pressure  contours  indicate  the  presence  of  a 
shock  between  the  two  airfoils,  with  a  Mach  number  ahead  of  the  shock  on  the  lower  air¬ 
foil  of  approximately  1.25.  The  costate  variables  shown  in  the  accompanying  figure,  on 
the  other  hand,  exhibit  a  shocklike  structure  in  a  location  that  corresponds  to  the  sonic  line 
in  the  flow  field.  However,  in  designing  for  other  objective  functions,  the  contours  of  the 
costate  variables  change  and  do  not  necessarily  show  such  a  clear  correspondence  with  the 
flow  field.  For  example,  if  the  cost  function  is  zero  at  the  design  point  in  an  unconstrained 
optimization,  the  costate  variables  are  all  zero  independent  of  the  flow  field,  due  to  the  ho¬ 
mogeneity  of  the  boundary  conditions. 
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The  initial  and  final  pressure  distributions  on  the  surface  of  the  airfoils  are  shown  in  Fig. 
19;  the  pressure  contours  after  three  design  cycles  are  shown  in  Fig.  20.  As  seen  from  Fig. 
19,  the  pressure  distribution  obtained  after  three  design  iterations  agrees  closely  with  that 
desired.  The  cost  function  has  been  reduced  over  3  orders  of  magnitude,  and  the  rms  of  the 
sensitivity  derivatives  has  also  been  reduced  over  3  orders  of  magnitude  after  the  second 
design  cycle.  The  final  pressure  distribution  on  the  aft  element  does  not  exhibit  the  strong 
shock  that  is  initially  present. 


Figure  IS.Contours  of  pressure  coefficient  and  'F,  for  two-element  configuration. 


Figure  20.  Final  airfoils  and  pressure  coefficient  for  two-element  airfoil. 

Viscous  Design  Examples 

For  the  first  viscous  case,  the  objective  is  to  maximize  the  lift  of  an  isolated  airfoil  by 
modifying  the  shape,  with  the  angle  of  attack  held  constant.  An  initial  computation  has 
been  performed  for  an  airfoil  at  a  Reynolds  number  of  5000  and  an  angle  of  attack  of  2° . 
The  mesh  used  for  this  computation  has  6951  nodes  of  which  128  lie  on  the  surface  of  the 
airfoil.  The  airfoil  geometry  is  described  by  using  a  12th-order  Bezier  representation  sim¬ 
ilar  to  that  described  earlier,  except  that  several  of  the  control  points  have  been  modified 
so  that  the  airfoil  is  no  longer  symmetric  (see  Fig.  21).  For  this  case,  nine  design  variables 
have  been  used.  These  correspond  to  the  y-coordinates  of  the  control  points  away  from  the 
immediate  vicinity  of  the  trailing  edge.  The  initial  lift  coefficient  is  0.0950,  and  the  initial 
drag  coefficient  is  0.0545.  After  three  design  cycles,  the  lift  has  been  increased  to  0.257 1 
and  5  of  the  nine  design  variables  have  hit  their  imposed  side  constraints.  Although  no 
constraint  or  objective  was  placed  on  the  drag,  the  drag  coefficient  has  dropped  to  0.0509. 
Note  that  for  this  case  both  the  initial  and  final  configurations  have  a  small  separated  re¬ 
gion  that  extends  over  the  last  25  percent  of  the  airfoil.  Despite  the  presence  of  separation, 
a  steady  flow  field  is  obtained.  In  the  event  of  unsteady  separation,  the  adjoint  approach  as 
described  would  not  be  applicable  because  the  steady-state  residual  is  assumed  to  be  zero 
and  is  used  as  a  constraint  for  the  optimization. 
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Figure  21.  Viscous  design  for  maximizing  lift. 

For  the  final  case  shown  in  Fig.  22,  the  objective  is  to  match  a  desired  pressure  distribu¬ 
tion  that  has  been  obtained  from  a  previous  analysis  of  an  NACA  0012  airfoil.  The  initial 
airfoil  geometry  has  been  obtained  by  simply  displacing  several  of  the  B-spline  control 
points  that  define  the  original  airfoil.  The  Reynolds  number  is  5000,  based  on  the  chord  of 
the  airfoil,  and  the  angle  of  attack  is  held  fixed  at  2° .  The  mesh  used  for  this  is  similar  to 
that  used  in  the  previous  test  case  and  has  approximately  7000  nodes.  For  the  current  test, 
the  cost  function  has  been  reduced  by  4.5  orders  of  magnitude  after  three  design  cycles, 
and  the  gradients  have  been  reduced  by  1.5  orders  of  magnitude  after  the  second  design 
cycle.  As  seen  in  Fig.  22,  the  target  pressure  distribution  is  obtained,  and  the  final  airfoil 
shape  is  that  of  an  NACA  0012  airfoil. 
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Figure  22.  Initial  and  final  pressure  distribution  for  viscous  flow. 

Discussion  and  Conclusions 

The  purpose  of  the  present  investigation  has  been  to  develop  and  analyze  the  continuous 
adjoint  approach  for  obtaining  sensitivity  derivatives  on  unstructured  grids  for  the  Euler 
and  Navier-Stokes  equations.  During  the  course  of  the  study,  several  drawbacks  have  been 
uncovered.  The  most  significant  is  the  need  for  accurate  second  derivatives  of  the  veloci¬ 
ties  required  for  computing  the  shape  sensitivity  derivatives  for  viscous  flows.  In  general, 
consistent  second  derivatives  cannot  be  obtained  with  spatially  second-order  accurate 
schemes.  This  problem  can  be  circumvented  by  mapping  the  domain  to  a  fixed  computa¬ 
tional  coordinate  system  as  is  usually  employed  for  structured  grids.  This  approach,  how¬ 
ever,  is  restrictive  in  its  generality  and  is  at  odds  with  the  flexibility  offered  by  unstruc¬ 
tured  grids.  The  absence  of  a  mapping  is  a  fundamental  difference  between  structured  and 
unstructured  grids.  The  requirement  for  second  derivatives  can  also  be  overcome  by  con¬ 
sidering  a  higher  order  discretization  of  the  flow  field,  so  that  consistent  second  deriva¬ 
tives  can  be  obtained.  However,  this  represents  a  significant  level  of  effort  because  the  en¬ 
tire  flow  field  must  be  computed  to  higher  order  accuracy.  It  appears  that  the  most 
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expedient  and  cost-effective  means  for  alleviating  this  problem  is  to  essentially  abandon 
the  purely  continuous  adjoint  approach  in  favor  of  a  more  discrete  approach,  as  described 
in  the  present  paper.  This  approach  has  the  added  benefit  that  the  contributions  to  the  sen¬ 
sitivity  derivatives  due  to  the  grid  may  be  included.  These  terms  do  not  naturally  appear  in 
the  continuous  framework  unless  a  mapping  to  a  fixed  computational  domain  is  first  em¬ 
ployed.  However,  it  is  shown  in  this  paper  that  these,  terms  are  critical  in  obtaining  accu¬ 
rate  derivatives  for  geometries  with  singularities. 

Another  possible  drawback  of  the  continuous  approach  is  that  a  restriction  is  placed  on 
the  allowable  cost  functions  that  can  be  used.  This  restriction  stems  from  the  need  for  a 
suitable  balance  between  terms  in  the  cost  function  and  corresponding  terms  from  the  re¬ 
sidual  that  are  used  to  eliminate  variations  in  pressure  or  velocity  gradients.  For  inviscid 
flows,  the  allowable  cost  functions  are  those  that  involve  only  the  pressure.  For  viscous 
flows,  an  additional  requirement  is  that  terms  from  the  entire  stress  tensor,  including  both 
the  prksure  and  viscous  terms,  must  be  included.  Although  this  limitation  does  not  appear 
to  occur  in  a  discrete  adjoint  approach,  the  full  implications  remain  unclear. 

The  continuous  adjoint  approach  requires  more  “up  front”  derivations  than  the  discrete 
approach  before  a  computer  implementation  can  be  pursued.  Also,  each  new  cost  function 
requires  a  certain  level  of  effort  to  not  only  arrive  at  the  appropriate  boundary  conditions 
but  to  determine  whether  boundary  conditions  can  even  be  obtained.  On  the  other  hand, 
the  continuous  approach  may  provide  insight  into  which  cost  functions  are  controllable. 
For  example,  in  the  case  of  inviscid  compressible  flow,  pressure  is  the  only  surviving  term 
in  the  boundary  flux  upon  application  of  the  flow  tangency  condition.  Therefore,  it  stands 
to  reason  that  only  cost  functions  that  involve  pressure  can  be  controlled.  In  the  discrete 
adjoint  approach,  new  cost  functions  are  more  easily  added  because  they  enter  the  prob¬ 
lem  only  through  the  right-hand  side  of  a  linear  system  of  equations.  After  a  subroutine 
has  been  written  to  evaluate  a  cost  function,  it  is  usually  a  simple  matter  to  obtain  all  the 
necessary  derivatives  by  differentiating  the  code  directly  using  the  chain  rule.  Further¬ 
more,  this  procedure  does  not  require  detailed  knowledge  of  the  equations  and  can  be  ac¬ 
complished  by  using  a  computational  tool  such  as  ADIFOR.® 

A  technique  is  presented  in  this  paper  that  is  derived  from  a  continuous  adjoint  approach 
but  appeals  to  the  discrete  approach  where  expedient.  A  discretization  of  the  adjoint  equa¬ 
tions  for  viscous  and  inviscid  flow  is  presented  that  corresponds  exactly  to  a  discrete  ad¬ 
joint  formulation  for  first-order  spatial  accuracy.  The  discretization  differs  from  the  dis¬ 
crete  adjoint  approach  for  higher  order  schemes  only  in  the  artificial  dissipation  terms. 
This  approach  is  simple  to  implement  and  yields  derivatives  that  are  reasonably  accurate 
in  comparison  with  finite-difference  calculations,  even  on  coarse  grids.  Alternatively,  the 
same  scheme  could  be  obtained  from  a  discrete  adjoint  point  of  view  by  appealing  to  the 
continuous  approach  for  making  suitable  approximations.  The  adjoint  approach  is  coupled 
with  an  optimization  algorithm  and  is  augmented  with  a  mesh  movement  strategy  for  re¬ 
structuring  the  mesh  in  response  to  surface  displacements.  The  mesh  movement  technique 
is  applicable  for  meshes  used  in  inviscid  computations  as  well  as  for  meshes  with  high  as¬ 
pect  ratio  triangles  typically  used  in  viscous  computations.  The  resulting  approach  has 
been  used  in  several  design  examples. 
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